Alejandro Naranjo Caraza
Orlando Almazán
import random
import plotly.graph_objects as go
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg, optimize
from AW_loadProblem import loadProblem
from mActiveSet import *
Imported: loadProblem (by AW).
I. Problema general de programación cuadrática¶
Un problema general de programación cuadrática (QP, por sus siglas en inglés) en su forma estándar de minimización puede formularse como:
$$ \begin{aligned} \min_{x \in \mathbb{R}^n} \quad & \frac{1}{2} x^\top G x + c^\top x \\ \text{sujeto a} \quad & A_i x = b_i, \quad \forall i \in \varepsilon, \\ & A_i x \leq b_i, \quad \forall i \in \mathcal{I}, \end{aligned} $$donde:
- $x \in \mathbb{R}^n$ es el vector de variables de decisión.
- $G \in \mathbb{R}^{n \times n}$ es una matriz simétrica definida positiva (o semidefinida positiva), que define el término cuadrático de la función objetivo.
- $c \in \mathbb{R}^n$ es el vector de coeficientes lineales en la función objetivo.
- $A \in \mathbb{R}^{m \times n}$ es la matriz de coeficientes de las restricciones.
- $b \in \mathbb{R}^m$ es el vector del lado derecho de las restricciones.
- $\varepsilon \subseteq \{1, \dots, m_1\}$ es el conjunto de índices correspondientes a restricciones de igualdad.
- $\mathcal{I} \subseteq \{1, \dots, m_2\}$ es el conjunto de índices correspondientes a restricciones de desigualdad.
Se asume que $\varepsilon \cap \mathcal{I} = \emptyset$ y que $\varepsilon \cup \mathcal{I} = \{1, \dots, m\}$. Es decir, cada restricción está clasificada de manera única como igualdad o desigualdad.
II. Implementación del problema¶
Definimos la clase QuadraticProgramming para representar un problema de programación cuadrática y evaluarlo, verificar su convexidad estricta y generar un punto factible.
__init__(...)¶
Inicializa el problema con:
G_matrix: matriz cuadrática $G$.c_vec: vector lineal $c$.A_matrix,b_vec: matriz y vector de restricciones (opcionales).eps_dic,l_dic: índices de las restricciones de igualdad ($\varepsilon$) y desigualdad ($\mathcal{I}$).
evaluate(x)¶
Calcula el valor de la función objetivo en un punto $x$ dado:
$$ f(x) = \frac{1}{2} x^\top G x + c^\top x $$generate_feasible()¶
Busca un punto factible resolviendo un problema lineal auxiliar con scipy.optimize.linprog.
- Si hay restricciones y están clasificados los índices (
epsyl), separa las filas correspondientes. - Si no hay restricciones, devuelve el vector cero.
- Si no encuentra solución, imprime un mensaje.
is_strict_convex()¶
Devuelve True si todos los autovalores de $G$ son positivos (es decir, $G$ es definida positiva), lo que implica que la función objetivo es estrictamente convexa:
class QuadraticProgramming:
def __init__(
self, G_matrix, c_vec, A_matrix=None, b_vec=None, eps_dic=None, l_dic=None
):
self.G = G_matrix
self.A = A_matrix
self.c = c_vec
self.b = b_vec
self.restrictions_eps = eps_dic
self.restrictions_l = l_dic
def evaluate(self, x):
res = (1 / 2) * x @ (self.G @ x) + self.c @ x
return res
def generate_feasible(self):
if self.A is not None and self.b is not None:
if self.restrictions_l or self.restrictions_eps:
A_eps = self.A[self.restrictions_eps, :]
A_l = self.A[self.restrictions_l, :]
b_eps = self.b[self.restrictions_eps]
b_l = self.b[self.restrictions_l]
c = np.zeros(A_eps.shape[1])
optim = optimize.linprog(c, A_ub=A_l, A_eq=A_eps, b_ub=b_l, b_eq=b_eps)
else:
optim = optimize.linprog(self.c, A_eq=self.A, b_eq=self.b)
if optim.success:
res = optim.x
else:
print("No feasible point found.")
print(optim.message)
res = None
else:
res = np.zeros(self.G.shape[1])
return res
def is_strict_convex(self):
return bool(np.all(np.linalg.eigvals(self.G) > 0))
Función schur_complement(problem)¶
Esta función resuelve un problema cuadrático con restricciones de igualdad usando el método del complemento de Schur.
Pasos clave:¶
- Extrae la submatriz de restricciones de igualdad $A$ y el vector $b$ desde el objeto
problem. - Calcula la descomposición de Cholesky de $G$:
$$ G = L^\top L $$ - Calcula $G^{-1}$ eficientemente usando $L$: $$ G^{-1} = (L^\top)^{-1} L^{-1} $$
- Define el sistema reducido: $$ M = L^{-1} A^\top $$
- Resuelve primero por los multiplicadores de Lagrange $\lambda$: $$ \lambda = -(M^\top M)^{-1} (b + A G^{-1} c) $$
- Luego obtiene $x$ óptimo: $$ x = -G^{-1} (c + A^\top \lambda) $$
Devuelve:¶
Una tupla (x, λ) con el punto óptimo y los multiplicadores.
Implementación¶
def schur_complement(problem):
if problem.A is None or problem.b is None:
raise ValueError("Schur complement expects constraints.")
if problem.restrictions_eps or problem.restrictions_l:
A = problem.A[problem.restrictions_eps, :]
b = problem.b[problem.restrictions_eps]
else:
A = problem.A
b = problem.b
A_transp = np.transpose(A)
G = problem.G
c = problem.c
L = np.transpose(linalg.cholesky(G))
L_transp = np.transpose(L)
L_inv = linalg.inv(L)
M = linalg.solve(L, A_transp)
M_transp = np.transpose(M)
G_inv = linalg.inv(L_transp) @ L_inv
lambd = linalg.solve(M_transp @ M, -b - A @ G_inv @ c)
x = -G_inv @ (c + A_transp @ lambd)
res = x, lambd
return res
Función null_space(problem)¶
Resuelve el problema usando el método del espacio nulo, útil cuando las restricciones de igualdad son de rango completo.
Pasos clave:¶
Verifica que $A$ tenga más columnas que filas y que tenga rango completo.
Aplica descomposición QR a $A^\top$: $$ A^\top = Q R $$ donde:
- $Y$ es la base del espacio columna de $A^\top$
- $Z$ es la base del espacio nulo de $A$
Descompone la solución en dos partes:
- $x_p$: solución particular que satisface $Ax = b$
- $x_h$: solución homogénea que vive en el espacio nulo de $A$
Calcula: $$ x_p = Y R^{-\top} b \\ x_h = Z (Z^\top G Z)^{-1} Z^\top (-c - G x_p) \\ x = x_p + x_h $$
Calcula los multiplicadores: $$ \lambda = R^{-1} (-Y^\top (c + G x)) $$
Devuelve:¶
Una tupla (x, λ) con la solución óptima y los multiplicadores.
Not: Ambas funciones suponen que el problema tiene solo restricciones de igualdad y que $G$ es simétrica definida positiva.
Implementación¶
def null_space(problem):
if problem.A is None or problem.b is None:
raise ValueError("Null space expects constraints.")
if problem.restrictions_eps or problem.restrictions_l:
A = problem.A[problem.restrictions_eps, :]
b = problem.b[problem.restrictions_eps]
else:
A = problem.A
b = problem.b
if A.shape[0] > A.shape[1]:
raise ValueError("Null space error. Matrix shape.")
rank_A = np.linalg.matrix_rank(A)
if rank_A < A.shape[0]:
raise ValueError("Null space error. Matrix not full row rank")
A_transp = np.transpose(A)
G = problem.G
c = problem.c
Q, R = linalg.qr(A_transp, mode="full")
Y = Q[:, :rank_A]
Z = Q[:, rank_A:]
R = R[:rank_A, :]
Z_transp = np.transpose(Z)
R_transp = np.transpose(R)
Y_transp = np.transpose(Y)
cp = linalg.solve(R_transp, b)
xp = Y @ cp
ch = linalg.solve(Z_transp @ G @ Z, -Z_transp @ (c + G @ xp))
x = xp + Z @ ch
lambd = linalg.solve(R, -Y_transp @ (c + G @ x))
res = x, lambd
return res
Método del Conjunto Activo¶
Entradas:¶
problem: instancia deQuadraticProgrammingx0: punto inicialW0: conjunto activo inicial (índices de restricciones activas)tol: tolerancia de convergenciamax_iter: número máximo de iteracionestrajectory: (opcional) si esTrue, guarda trayectoria de $x_k$ y $q_k$
Salidas:¶
- Solución óptima
x* - Conjunto activo final
W* - (Opcional) trayectoria de iteraciones
Estructura:¶
active_set_branch_1: busca dirección $d_k$ resolviendo un subproblema QP reducido:- Usa Schur complement o espacio nulo según el caso.
- Calcula paso $\alpha$ y actualiza $x_{k+1} = x_k + \alpha d_k$
- Si alguna restricción desigual se activa, se añade a $W_k$.
active_set_branch_2: revisa los multiplicadores de Lagrange $\lambda_k$:- Si alguno $\lambda_j < 0$ y $j \notin \varepsilon$, elimina $j$ de $W_k$
- Si todos $\lambda_j \geq 0$, se ha encontrado solución.
Cicla entre rama 1 y rama 2 hasta converger o alcanzar max_iter.
Implementación¶
def active_set_branch_1(
problem,
x0,
W0,
zeros,
x_trajectory=None,
q_trajectory=None,
tol=10e-10,
max_iter=100,
):
k = 0
xk = x0
Wk = W0
next_iter = True
if x_trajectory is not None:
x_trajectory.append(xk)
if q_trajectory is not None:
q_trajectory.append(problem.evaluate(xk))
while next_iter and k < max_iter:
print(f"\n\033[4mBranch 1; k = {k}\033[0m")
print(f"W{k} = {Wk}")
gk = problem.G @ xk + problem.c
if Wk:
Ak = problem.A[Wk, :]
zerosk = zeros[Wk]
problem_k = QuadraticProgramming(problem.G, gk, Ak, zerosk)
schur_condition_1 = problem_k.is_strict_convex()
schur_condition_2 = (
abs(problem_k.A.shape[1] - problem_k.A.shape[0])
> problem_k.A.shape[1] / 2
)
if schur_condition_1 and schur_condition_2:
dk, lambdk = schur_complement(problem_k)
print("schur complement used")
else:
dk, lambdk = null_space(problem_k)
print("null space used")
else:
problem_k = QuadraticProgramming(problem.G, gk)
dk, lambdk = non_restriction_optimization(problem_k)
print("non restriction optimization")
dk_norm = np.linalg.norm(dk, ord=np.inf)
print(f"norm(d{k}) = {np.round(dk_norm, decimals=3)}")
if dk_norm > tol:
Wk_set = set(Wk)
temp_set_l = [i for i in problem.restrictions_l if i not in Wk_set]
Al = problem.A[temp_set_l, :]
bl = problem.b[temp_set_l]
alpha_condition = (Al @ dk > 0) & (bl - Al @ xk > 0)
temp_set_l = np.array(temp_set_l)[alpha_condition].tolist()
if not temp_set_l:
alpha = 1
print("No alpha tilde.")
else:
Al = Al[alpha_condition, :]
bl = np.extract(alpha_condition, bl)
alpha_vec = (bl - Al @ xk) / (Al @ dk)
index_min = np.argmin(alpha_vec)
alpha_min = alpha_vec[index_min]
j_min = int(temp_set_l[index_min])
print(
f"alpha tilde = {
np.round(alpha_min, decimals=8)}"
)
alpha = min(alpha_min, 1.0)
xk = xk + alpha * dk
if x_trajectory is not None:
x_trajectory.append(xk)
if q_trajectory is not None:
q_trajectory.append(problem.evaluate(xk))
print(f"x{k+1} = {np.round(xk, decimals=3)}")
print(f"q{k+1} = {np.round(problem.evaluate(xk), decimals=3)}")
if alpha != 1:
print(f"W{k+1} = W{k} U [{j_min}]")
Wk.append(j_min)
k += 1
else:
next_iter = False
print(f"W{k+1} = W{k}")
else:
next_iter = False
if next_iter and k == max_iter:
raise MaxIterationsExceeded(
f"Branch 1 maximum number of iterations ({max_iter}) exceeded."
)
else:
res = xk, Wk, lambdk
return res
def active_set_branch_2(problem, xk, Wk, lambdk):
print("\n\033[4mBranch 2\033[0m")
min_mu = 0
min_j = None
for i in range(len(Wk)):
if Wk[i] not in problem.restrictions_eps and lambdk[i] < min_mu:
min_mu = lambdk[i]
min_j = Wk[i]
if min_mu < 0:
Wk.remove(min_j)
next_iter = True
print(f"m{min_j} = {min_mu} < 0. Back to Branch 1.")
else:
next_iter = False
print("mj >= 0 forall j. Solution found.")
return Wk, next_iter
def active_set(problem, x0, W0, tol=10e-9, max_iter=100, trajectory=False):
if not isinstance(problem, QuadraticProgramming):
raise TypeError("Expected Quadratic Programming Instance.")
if not linalg.issymmetric(problem.G):
raise ValueError("Optimization matrix is not symetric")
if not set(problem.restrictions_eps).issubset(W0):
raise ValueError("Error. W0 not subset of equality restrictions")
if trajectory:
x_trajectory = []
q_trajectory = []
else:
x_trajectory = None
q_trajectory = None
xk = np.copy(x0)
Wk = W0.copy()
zeros = np.zeros(problem.A.shape[0])
iter = 0
next_iter = True
while next_iter and iter < max_iter:
print("\n-----------------------------------------------------------")
print(f"ITERATION {iter}")
print("-----------------------------------------------------------")
xk, Wk, lambdk = active_set_branch_1(
problem, xk, Wk, zeros, x_trajectory, q_trajectory
)
Wk, next_iter = active_set_branch_2(problem, xk, Wk, lambdk)
iter += 1
print("\n-----------------------------------------------------------")
if next_iter and iter == max_iter:
res = None
print("Outer layer maximum iterations reached. No solution found.")
print("-----------------------------------------------------------\n")
else:
if trajectory:
x_trajectory = np.array(x_trajectory)
q_trajectory = np.array(q_trajectory)
res = xk, Wk, x_trajectory, q_trajectory
else:
res = xk, Wk
print("ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND.")
print("-----------------------------------------------------------\n")
print(f"x* = {np.round(xk, decimals=3)}\nW* = {Wk}")
return res
III. Problema 1¶
En el primer problema se pide minimizar la función objetivo:
$$ q(\vec{x}) = (x - 1)^2 + (y - 2.5)^2 $$sujeta a las siguientes restricciones:
$$ \begin{aligned} -x + 2y - 2 &\leq 0 \\ x + 2y - 6 &\leq 0 \\ x - 2y - 2 &\leq 0 \\ -x &\leq 0 \\ -y &\leq 0 \end{aligned} $$Primera prueba de optimización¶
Partimos del siguiente punto inicial:
$$ \vec{x}_0 = \begin{pmatrix} 2 \\ 0 \end{pmatrix} $$y del conjunto activo inicial:
$$ W_0 = \{3\} \subset \{1, 2, 3, 4, 5\} $$Nota: En la implementación del algoritmo, las restricciones están indexadas desde cero (0, 1, 2, ...), por lo que al ejecutar la función se toma $W_0 = \{2\}$.
G_mat = np.array([[2, 0], [0, 2]])
A_mat = np.array([[-1, 2], [1, 2], [1, -2], [-1, 0], [0, -1]])
c_vec = np.array([-2, -5])
b_vec = np.array([2, 6, 2, 0, 0])
l_dic = [0, 1, 2, 3, 4]
eps_dic = []
problem = QuadraticProgramming(G_mat, c_vec, A_mat, b_vec, eps_dic, l_dic)
W_0 = [2]
x0 = np.array([2, 0])
res = active_set(problem, x0, W_0, trajectory=True, tol=1e-9)
xf, Wf, x_traj, q_traj = res
----------------------------------------------------------- ITERATION 0 ----------------------------------------------------------- Branch 1; k = 0 W0 = [2] null space used norm(d0) = 0.2 alpha tilde = 10.0 x1 = [2.2 0.1] q1 = -0.05 W1 = W0 Branch 2 m2 = -2.3999999999999995 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 1 ----------------------------------------------------------- Branch 1; k = 0 W0 = [] non restriction optimization norm(d0) = 2.4 alpha tilde = 0.66666667 x1 = [1.4 1.7] q1 = -6.45 W1 = W0 U [0] Branch 1; k = 1 W1 = [0] null space used norm(d1) = 0.0 Branch 2 mj >= 0 forall j. Solution found. ----------------------------------------------------------- ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND. ----------------------------------------------------------- x* = [1.4 1.7] W* = [0]
Resultados¶
print("Optimum reached at x* =",xf)
print("Function at x* is q(x*) =",problem.evaluate(xf))
print("Final active set is W* =",Wf)
Optimum reached at x* = [1.4 1.7] Function at x* is q(x*) = -6.45 Final active set is W* = [0]
x = np.linspace(0, 6, 100)
y = np.linspace(0, 6, 100)
x, y = np.meshgrid(x, y)
z = 0.5 * (problem.G[0,0]*x**2 + 2*problem.G[0,1]*x*y + problem.G[1,1]*y**2) + problem.c[0]*x + problem.c[1]*y
points_x = x_traj[:,0]
points_y = x_traj[:,1]
points_z = q_traj
fig = go.Figure()
fig.add_trace(go.Surface(z=z, x=x, y=y, colorscale='Viridis', opacity=0.9))
fig.add_trace(go.Scatter3d(
x=points_x,
y=points_y,
z=points_z,
mode='markers+text',
marker=dict(size=6, color='red'),
text=[f"P{i}" for i in range (points_x.size)],
textposition="top center",
textfont=dict(color="red")
))
fig.update_layout(
title='Trajectory plot',
scene=dict(
xaxis_title='X',
yaxis_title='Y',
zaxis_title='Z'
),
margin=dict(l=0, r=0, b=0, t=30)
)
fig.show()
A = problem.A
b = problem.b
xtraj = x_traj[:, 0]
ytraj = x_traj[:, 1]
x = np.linspace(0, 6, 900)
y = np.linspace(0, 6, 900)
X, Y = np.meshgrid(x, y)
points = np.vstack((X.ravel(), Y.ravel())).T
feasible = np.all(A @ points.T <= b[:, np.newaxis], axis=0)
X_feasible = X.ravel()[feasible]
Y_feasible = Y.ravel()[feasible]
plt.figure(figsize=(8, 8))
plt.scatter(X_feasible, Y_feasible, color='grey', s=0.5, alpha=0.075)
plt.scatter(xtraj, ytraj, color='red', s=30)
for i in range (len(xtraj)):
plt.text(xtraj[i], ytraj[i], f"P{i}", fontsize=12, ha='right')
plt.xlim(0, 5)
plt.ylim(0, 2.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Feasible Region and Point Trajectory')
plt.grid(True)
plt.show()
Usamos la biblioteca CVXPY para comparar nuestros resultados.
# Compare the CVXPY solution.
n = int(problem.G.shape[1])
x = cp.Variable(n)
objective = cp.Minimize((1 / 2) * cp.quad_form(x, problem.G) + problem.c.T @ x)
constraints = [problem.A[problem.restrictions_l, :] @ x <= problem.b[problem.restrictions_l]]
prob = cp.Problem(objective, constraints)
prob.solve()
solution = x.value
print("Method used:",prob.solver_stats.solver_name)
print("Our optimum value:", round(float(problem.evaluate(xf)),2))
print("Optimum value with cvxpy:",round(problem.evaluate(x.value),2))
print("Solution x with active set method:",np.round(xf,2))
print("Solution x with cvxpy: ",np.round(solution,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-x.value)/(np.linalg.norm(xf)+np.linalg.norm(x.value)),2))
Method used: OSQP Our optimum value: -6.45 Optimum value with cvxpy: -6.45 Solution x with active set method: [1.4 1.7] Solution x with cvxpy: [1.4 1.7] Normalized euclidean distance: 0.0
Observaciones¶
- Las soluciones obtenidas con CVXPY y nuestro implementación del Método del Conjunto Activo coinciden.
- El Método de Conjunto Activo encuentra nuestra solución en dos iteraciones (trés puntos en total).
- CVXPY utiliza el método OSQP, basado en separación de operadores utilizando ADMM (Alternating Direction Method of Multipliers).
Segunda prueba de optimización¶
Partimos del siguiente punto inicial:
$$ \vec{x}_0 = \begin{pmatrix} 3 \\ 1 \end{pmatrix} $$y del conjunto activo inicial:
$$ W_0 = \emptyset\subset \{1, 2, 3, 4, 5\} $$G_mat = np.array([[2, 0], [0, 2]])
A_mat = np.array([[-1, 2], [1, 2], [1, -2], [-1, 0], [0, -1]])
c_vec = np.array([-2, -5])
b_vec = np.array([2, 6, 2, 0, 0])
l_dic = [0, 1, 2, 3, 4]
eps_dic = []
problem = QuadraticProgramming(G_mat, c_vec, A_mat, b_vec, eps_dic, l_dic)
W_0 = [1]
x0 = np.array([3, 1])
res = active_set(problem, x0, W_0, trajectory=True, tol=1e-9)
xf, Wf, x_traj, q_traj = res
----------------------------------------------------------- ITERATION 0 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1] null space used norm(d0) = 2.2 alpha tilde = 0.68181818 x1 = [1.5 1.75] q1 = -6.437 W1 = W0 U [0] Branch 1; k = 1 W1 = [1, 0] null space used norm(d1) = 0.0 Branch 2 m1 = -0.1250000000000004 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 1 ----------------------------------------------------------- Branch 1; k = 0 W0 = [0] null space used norm(d0) = 0.1 alpha tilde = 15.0 x1 = [1.4 1.7] q1 = -6.45 W1 = W0 Branch 2 mj >= 0 forall j. Solution found. ----------------------------------------------------------- ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND. ----------------------------------------------------------- x* = [1.4 1.7] W* = [0]
Resultados¶
print("Optimum reached at x* =",xf)
print("Function at x* is q(x*) =",round(problem.evaluate(xf),2))
print("Final active set is W* =",Wf)
Optimum reached at x* = [1.4 1.7] Function at x* is q(x*) = -6.45 Final active set is W* = [0]
x = np.linspace(0, 6, 100)
y = np.linspace(0, 6, 100)
x, y = np.meshgrid(x, y)
z = 0.5 * (problem.G[0,0]*x**2 + 2*problem.G[0,1]*x*y + problem.G[1,1]*y**2) + problem.c[0]*x + problem.c[1]*y
points_x = x_traj[:,0]
points_y = x_traj[:,1]
points_z = q_traj
fig = go.Figure()
fig.add_trace(go.Surface(z=z, x=x, y=y, colorscale='Viridis', opacity=0.9))
fig.add_trace(go.Scatter3d(
x=points_x,
y=points_y,
z=points_z,
mode='markers+text',
marker=dict(size=6, color='red'),
text=[f"P{i}" for i in range (len(points_x))],
textposition="top center",
textfont=dict(color="red")
))
fig.update_layout(
title='Trajectory plot',
scene=dict(
xaxis_title='X',
yaxis_title='Y',
zaxis_title='Z'
),
margin=dict(l=0, r=0, b=0, t=30)
)
fig.show()
A = problem.A
b = problem.b
xtraj = x_traj[:, 0]
ytraj = x_traj[:, 1]
x = np.linspace(0, 6, 900)
y = np.linspace(0, 6, 900)
X, Y = np.meshgrid(x, y)
points = np.vstack((X.ravel(), Y.ravel())).T
feasible = np.all(A @ points.T <= b[:, np.newaxis], axis=0)
X_feasible = X.ravel()[feasible]
Y_feasible = Y.ravel()[feasible]
plt.figure(figsize=(8, 8))
plt.scatter(X_feasible, Y_feasible, color='grey', s=0.5, alpha=0.075)
plt.scatter(xtraj, ytraj, color='red', s=30)
for i in range (len(xtraj)):
plt.text(xtraj[i], ytraj[i], f"P{i}", fontsize=12, ha='right')
plt.xlim(0, 5)
plt.ylim(0, 2.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Feasible Region and Point Trajectory')
plt.grid(True)
plt.show()
# Compare with CVXPY solution
n = int(problem.G.shape[1])
x = cp.Variable(n)
objective = cp.Minimize((1 / 2) * cp.quad_form(x, problem.G) + problem.c.T @ x)
constraints = [problem.A[problem.restrictions_l, :] @ x <= problem.b[problem.restrictions_l]]
prob = cp.Problem(objective, constraints)
prob.solve()
solution = x.value
print("Our optimum value:", round(float(problem.evaluate(xf)),2))
print("Optimum value obtained with cvxpy:",round(problem.evaluate(x.value),2))
print("Solution x with active set method:",np.round(xf,2))
print("Solution x with cvxpy: ",np.round(solution,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-x.value)/(np.linalg.norm(xf)+np.linalg.norm(x.value)),2))
Our optimum value: -6.45 Optimum value obtained with cvxpy: -6.45 Solution x with active set method: [1.4 1.7] Solution x with cvxpy: [1.4 1.7] Normalized euclidean distance: 0.0
Observaciones¶
- Obtenemos una solución idéntica
- En este caso, comenzamos en el interior de nuestra región factible. El algoritmo se mueve a la frontera de la región y finalmente optimiza en la frontera.
IV. Problema - Klee Minty¶
Para este problema, se considera la matriz $ G = 10^{-4} I $ y el problema cuadrático de Klee-Minty definido por:
$$ \min_{\vec{x} \in \mathbb{R}^n} \left\{ \frac{1}{2} \vec{x}^\top G \vec{x} - \sum_{i=1}^n x_i \right\} $$sujeto a las siguientes restricciones:
$$ \begin{aligned} x_1 &\leq 1 \\ 2 \sum_{j=1}^{i-1} x_j + x_i &\leq 2^{i-1}, \quad i = 2, \dots, n \\ x_1, \dots, x_n &\geq 0 \end{aligned} $$Sea $n = 20$. Este problema impone $n$ restricciones de no negatividad, una para cada variable. Para aplicar el método de optimización (por ejemplo, un método de conjunto activo o de proyecciones), se debe construir un punto inicial $\vec{x}_{0}$ que cumpla con un conjunto activo inicial $W_0$. Este conjunto debe ser un subconjunto aleatorio formado por 8 de las últimas 10 restricciones de no negatividad.
A partir del punto $\vec{x}_{0}$ y del conjunto activo $W_0$, se debe ejecutar el método seleccionado y determinar:
- El conjunto activo inicial $W_0$,
- El número de iteraciones necesarias,
- El valor óptimo alcanzado.
def klee_minty(n):
G = np.identity(n) * 1e-4
c = np.ones(n) * -1
b = 2 ** np.arange(1, n + 1) - 1
b = np.concatenate((b, np.zeros(n)))
L = np.tril(np.ones((n, n))) + np.tril(np.ones((n, n)), k=-1)
Identity = np.identity(n)
A = np.vstack((L, -Identity))
restrictions_eps = []
restrictions_l = [i for i in range(n)]
problem = QuadraticProgramming(G, c, A, b, restrictions_eps, restrictions_l)
return problem
problem_km = klee_minty(20)
k = random.randint(1, 8)
W0 = (
np.random.choice(problem_km.restrictions_l[-8:], k, replace=False)
.astype(int)
.tolist()
)
problem_km_0=klee_minty(20)
problem_km_0.restrictions_eps = list(set(problem_km_0.restrictions_eps + W0))
problem_km_0.restrictions_l = [
i for i in problem_km_0.restrictions_l if i not in W0
]
x0 = problem_km_0.generate_feasible()
res = active_set(problem_km, x0, W0, max_iter=100000,
tol=1e-9, trajectory=True)
xf, Wf, x_traj, q_traj = res
----------------------------------------------------------- ITERATION 0 ----------------------------------------------------------- Branch 1; k = 0 W0 = [18, 12, 17, 14, 13, 16, 19, 15] schur complement used norm(d0) = 23892.447 No alpha tilde. x1 = [ -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 -952.852 31059.447 -22867.447 39251.447 -6483.447 72019.447 59052.553 203091.447 321196.553] q1 = 7123543.321 W1 = W0 Branch 2 m19 = -31.119655324675243 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 1 ----------------------------------------------------------- Branch 1; k = 0 W0 = [18, 12, 17, 14, 13, 16, 15] schur complement used norm(d0) = 311196.553 alpha tilde = 0.51647119 x1 = [1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 8.167000e+03 2.500000e+01 1.635900e+04 1.640900e+04 4.912700e+04 8.194500e+04 1.801990e+05 1.604725e+05] q1 = 2885040.535 W1 = W0 U [1] Branch 1; k = 1 W1 = [18, 12, 17, 14, 13, 16, 15, 1] schur complement used norm(d1) = 150472.5 alpha tilde = 0.00136048 x2 = [7.14000000e-01 1.57100000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00 8.13785700e+03 5.41430000e+01 1.63298570e+04 1.64381430e+04 4.90978570e+04 8.19741430e+04 1.80169857e+05 1.60267786e+05] q2 = 2881524.022 W2 = W1 U [2] Branch 1; k = 2 W2 = [18, 12, 17, 14, 13, 16, 15, 1, 2] schur complement used norm(d2) = 150267.786 alpha tilde = 0.00219133 x3 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 8.091000e+03 1.010000e+02 1.628300e+04 1.648500e+04 4.905100e+04 8.202100e+04 1.801230e+05 1.599385e+05] q3 = 2875877.999 W3 = W2 U [3] Branch 1; k = 3 W3 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3] null space used norm(d3) = 149938.5 alpha tilde = 0.00393895 x4 = [6.000000e-01 1.800000e+00 2.200000e+00 5.800000e+00 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01 8.007000e+03 1.850000e+02 1.619900e+04 1.656900e+04 4.896700e+04 8.210500e+04 1.800390e+05 1.593479e+05] q4 = 2865783.069 W4 = W3 U [4] Branch 1; k = 4 W4 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4] null space used norm(d4) = 149347.9 alpha tilde = 0.00452233 x5 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 1.700000e+01 1.700000e+01 1.700000e+01 1.700000e+01 1.700000e+01 1.700000e+01 1.700000e+01 7.911000e+03 2.810000e+02 1.610300e+04 1.666500e+04 4.887100e+04 8.220100e+04 1.799430e+05 1.586725e+05] q5 = 2854288.422 W5 = W4 U [0] Branch 1; k = 5 W5 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0] null space used norm(d5) = 148672.5 alpha tilde = 0.00265012 x6 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01 7.855000e+03 3.370000e+02 1.604700e+04 1.672100e+04 4.881500e+04 8.225700e+04 1.798870e+05 1.582785e+05] q6 = 2847607.025 W6 = W5 U [5] Branch 1; k = 6 W6 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5] null space used norm(d6) = 148278.5 alpha tilde = 0.01253722 x7 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 4.300000e+01 4.300000e+01 4.300000e+01 4.300000e+01 4.300000e+01 7.591000e+03 6.010000e+02 1.578300e+04 1.698500e+04 4.855100e+04 8.252100e+04 1.796230e+05 1.564195e+05] q7 = 2816324.832 W7 = W6 U [6] Branch 1; k = 7 W7 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6] null space used norm(d7) = 146419.5 alpha tilde = 0.02022272 x8 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 8.500000e+01 8.500000e+01 8.500000e+01 8.500000e+01 8.500000e+01 7.171000e+03 1.021000e+03 1.536300e+04 1.740500e+04 4.813100e+04 8.294100e+04 1.792030e+05 1.534585e+05] q8 = 2767320.956 W8 = W7 U [7] Branch 1; k = 8 W8 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7] null space used norm(d8) = 143458.5 alpha tilde = 0.03387042 x9 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 1.710000e+02 1.710000e+02 1.710000e+02 6.483000e+03 1.709000e+03 1.467500e+04 1.809300e+04 4.744300e+04 8.362900e+04 1.785150e+05 1.485995e+05] q9 = 2689092.17 W9 = W8 U [8] Branch 1; k = 9 W9 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8] null space used norm(d9) = 138599.5 alpha tilde = 0.05212862 x10 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02 3.410000e+02 3.410000e+02 5.463000e+03 2.729000e+03 1.365500e+04 1.911300e+04 4.642300e+04 8.464900e+04 1.774950e+05 1.413745e+05] q10 = 2577795.388 W10 = W9 U [9] Branch 1; k = 10 W10 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8, 9] null space used norm(d10) = 131374.5 alpha tilde = 0.07419248 x11 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02 6.830000e+02 6.830000e+02 4.095000e+03 4.097000e+03 1.228700e+04 2.048100e+04 4.505500e+04 8.601700e+04 1.761270e+05 1.316275e+05] q11 = 2437189.527 W11 = W10 U [10] Branch 1; k = 11 W11 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8, 9, 10] null space used norm(d11) = 121627.5 alpha tilde = 0.08130563 x12 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02 6.830000e+02 1.365000e+03 2.731000e+03 5.461000e+03 1.092300e+04 2.184500e+04 4.369100e+04 8.738100e+04 1.747630e+05 1.217385e+05] q12 = 2305886.147 W12 = W11 U [11] Branch 1; k = 12 W12 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8, 9, 10, 11] null space used norm(d12) = 111738.5 No alpha tilde. x13 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 4.36910e+04 8.73810e+04 1.74763e+05 1.00000e+04] q13 = 1681611.528 W13 = W12 Branch 2 m0 = -22.302900000007853 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 2 ----------------------------------------------------------- Branch 1; k = 0 W0 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] null space used norm(d0) = 6110.384 No alpha tilde. x1 = [ -3054.192 6111.384 -6107.384 6115.384 -6099.384 6131.384 -6067.384 6195.384 -5939.384 6451.384 -5427.384 7475.384 -3379.384 11571.384 4812.616 27955.384 37580.616 93491.384 168652.616 10000. ] q1 = 1647541.709 W1 = W0 Branch 2 m16 = -17.79030821918008 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 3 ----------------------------------------------------------- Branch 1; k = 0 W0 = [18, 12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] null space used norm(d0) = 40243.785 alpha tilde = 0.95655556 x1 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02 6.830000e+02 1.365000e+03 2.731000e+03 5.461000e+03 1.092300e+04 2.184500e+04 2.138810e+04 1.319868e+05 1.301572e+05 1.000000e+04] q1 = 1440535.133 W1 = W0 U [0] Branch 1; k = 1 W1 = [18, 12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0] null space used norm(d1) = 1717.533 No alpha tilde. x2 = [1.00000000e+00 1.00000000e+00 3.00000000e+00 5.00000000e+00 1.10000000e+01 2.10000000e+01 4.30000000e+01 8.50000000e+01 1.71000000e+02 3.41000000e+02 6.83000000e+02 1.36500000e+03 2.73100000e+03 5.46100000e+03 1.09230000e+04 2.18450000e+04 2.05293330e+04 1.33704333e+05 1.28439667e+05 1.00000000e+04] q2 = 1440203.266 W2 = W1 Branch 2 m18 = -11.843966666666557 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 4 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0] null space used norm(d0) = 118439.667 alpha tilde = 0.48889167 x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 4.36910e+04 8.73810e+04 7.05355e+04 1.00000e+04] q1 = 507496.557 W1 = W0 U [16] Branch 1; k = 1 W1 = [12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16] null space used norm(d1) = 60535.5 No alpha tilde. x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03 2.7310e+03 5.4610e+03 1.0923e+04 2.1845e+04 4.3691e+04 8.7381e+04 1.0000e+04 1.0000e+04] q2 = 324269.219 W2 = W1 Branch 2 m3 = -10.649699999991908 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 5 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 17, 14, 13, 15, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16] null space used norm(d0) = 3736.737 No alpha tilde. x1 = [ 1.0000000e+00 1.0000000e+00 3.0000000e+00 -1.8633680e+03 3.7477370e+03 -3.7157370e+03 3.7797370e+03 -3.6517370e+03 3.9077370e+03 -3.3957370e+03 4.4197370e+03 -2.3717370e+03 6.4677370e+03 1.7242630e+03 1.4659737e+04 1.8108263e+04 4.7427737e+04 8.3644263e+04 1.0000000e+04 1.0000000e+04] q1 = 314320.438 W1 = W0 Branch 2 m15 = -8.054131578943784 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 6 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 17, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16] null space used norm(d0) = 18331.437 alpha tilde = 0.95806252 x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 5.46100e+03 1.09230e+04 1.11953e+04 6.49904e+04 6.60816e+04 1.00000e+04 1.00000e+04] q1 = 269634.821 W1 = W0 U [3] Branch 1; k = 1 W1 = [12, 17, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16, 3] null space used norm(d1) = 750.6 No alpha tilde. x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03 2.7310e+03 5.4610e+03 1.0923e+04 1.0820e+04 6.5741e+04 6.5331e+04 1.0000e+04 1.0000e+04] q2 = 269571.438 W2 = W1 Branch 2 m17 = -5.533100000000015 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 7 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16, 3] null space used norm(d0) = 55331.0 alpha tilde = 0.49813848 x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 4.36910e+04 3.77685e+04 1.00000e+04 1.00000e+04] q1 = 63432.741 W1 = W0 U [15] Branch 1; k = 1 W1 = [12, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16, 3, 15] null space used norm(d1) = 27768.5 No alpha tilde. x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03 2.7310e+03 5.4610e+03 1.0923e+04 2.1845e+04 4.3691e+04 1.0000e+04 1.0000e+04 1.0000e+04] q2 = 24878.261 W2 = W1 Branch 2 m0 = -4.826500000005277 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 8 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 16, 3, 15] null space used norm(d0) = 1485.077 No alpha tilde. x1 = [ -741.538 1486.077 -1482.077 1490.077 -1474.077 1506.077 -1442.077 1570.077 -1314.077 1826.077 -802.077 2850.077 1245.923 6946.077 9437.923 23330.077 42205.923 10000. 10000. 10000. ] q1 = 23086.33 W1 = W0 Branch 2 m14 = -3.7189615384644465 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 9 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 16, 3, 15] null space used norm(d0) = 8434.758 alpha tilde = 0.96836482 x1 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03 2.7310e+03 5.4610e+03 6.0965e+03 3.1498e+04 3.4038e+04 1.0000e+04 1.0000e+04 1.0000e+04] q1 = 13827.628 W1 = W0 U [0] Branch 1; k = 1 W1 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 16, 3, 15, 0] null space used norm(d1) = 261.444 No alpha tilde. x2 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02 6.8300000e+02 1.3650000e+03 2.7310000e+03 5.4610000e+03 5.9657780e+03 3.1759444e+04 3.3776556e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q2 = 13819.938 W2 = W1 Branch 2 m16 = -2.3776555555555046 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 10 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 3, 15, 0] null space used norm(d0) = 23776.556 alpha tilde = 0.52123006 x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 2.13835e+04 1.00000e+04 1.00000e+04 1.00000e+04] q1 = -25396.709 W1 = W0 U [14] Branch 1; k = 1 W1 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 3, 15, 0, 14] null space used norm(d1) = 11383.5 No alpha tilde. x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03 2.7310e+03 5.4610e+03 1.0923e+04 2.1845e+04 1.0000e+04 1.0000e+04 1.0000e+04 1.0000e+04] q2 = -31875.913 W2 = W1 Branch 2 m3 = -1.9116999999972746 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 11 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 15, 0, 14] null space used norm(d0) = 780.286 No alpha tilde. x1 = [ 1.0000000e+00 1.0000000e+00 3.0000000e+00 -3.8514300e+02 7.9128600e+02 -7.5928600e+02 8.2328600e+02 -6.9528600e+02 9.5128600e+02 -4.3928600e+02 1.4632860e+03 5.8471400e+02 3.5112860e+03 4.6807140e+03 1.1703286e+04 2.1064714e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q1 = -32248.831 W1 = W0 Branch 2 m13 = -1.3403571428555356 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 12 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 15, 0, 14] null space used norm(d0) = 3065.92 alpha tilde = 0.99256168 x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 3.54930e+03 1.47464e+04 1.80216e+04 1.00000e+04 1.00000e+04 1.00000e+04 1.00000e+04] q1 = -33539.541 W1 = W0 U [3] Branch 1; k = 1 W1 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 15, 0, 14, 3] null space used norm(d1) = 22.156 No alpha tilde. x2 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02 6.8300000e+02 1.3650000e+03 2.7310000e+03 3.5382220e+03 1.4768556e+04 1.7999444e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q2 = -33539.596 W2 = W1 Branch 2 m15 = -0.7999444444444724 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 13 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 14, 3] null space used norm(d0) = 7999.444 alpha tilde = 0.60090979 x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03 2.73100e+03 5.46100e+03 1.09230e+04 1.31925e+04 1.00000e+04 1.00000e+04 1.00000e+04 1.00000e+04] q1 = -38381.511 W1 = W0 U [13] Branch 1; k = 1 W1 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 14, 3, 13] null space used norm(d1) = 3192.5 No alpha tilde. x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03 2.7310e+03 5.4610e+03 1.0923e+04 1.0000e+04 1.0000e+04 1.0000e+04 1.0000e+04 1.0000e+04] q2 = -38891.114 W2 = W1 Branch 2 m0 = -0.45730000000121807 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 14 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 14, 3, 13] null space used norm(d0) = 160.456 No alpha tilde. x1 = [-7.9228000e+01 1.6145600e+02 -1.5745600e+02 1.6545600e+02 -1.4945600e+02 1.8145600e+02 -1.1745600e+02 2.4545600e+02 1.0544000e+01 5.0145600e+02 5.2254400e+02 1.5254560e+03 2.5705440e+03 5.6214560e+03 1.0762544e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q1 = -38909.458 W1 = W0 Branch 2 m10 = -0.2894894736850741 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 15 ----------------------------------------------------------- Branch 1; k = 0 W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13] null space used norm(d0) = 350.08 alpha tilde = 0.9930721 x1 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02 4.2894400e+02 1.8731110e+03 2.2228890e+03 5.9691110e+03 1.0414889e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q1 = -38946.486 W1 = W0 U [0] Branch 1; k = 1 W1 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13, 0] null space used norm(d1) = 2.359 No alpha tilde. x2 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02 4.2776500e+02 1.8754710e+03 2.2205290e+03 5.9714710e+03 1.0412529e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q2 = -38946.487 W2 = W1 Branch 2 m12 = -0.11026470588274505 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 16 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13, 0] null space used norm(d0) = 416.093 alpha tilde = 1.22681782 x1 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02 6.3581100e+02 1.4593770e+03 2.2829430e+03 6.2627360e+03 1.0121264e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q1 = -38965.987 W1 = W0 Branch 2 m0 = -0.018681132075438582 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 17 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13] null space used norm(d0) = 9.876 alpha tilde = 8.16878223 x1 = [-3.9380000e+00 1.0876000e+01 -6.8760000e+00 1.4876000e+01 1.1240000e+00 3.0876000e+01 3.3124000e+01 9.4876000e+01 1.6112400e+02 3.5087600e+02 6.3171200e+02 1.4577000e+03 2.2836890e+03 6.2629220e+03 1.0121078e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04] q1 = -38966.033 W1 = W0 Branch 2 m14 = -0.012107780548628416 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 18 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 2, 4, 5, 6, 7, 8, 9, 11, 3, 13] null space used norm(d0) = 121.078 alpha tilde = 6.06203426 x1 = [-3.717000e+00 1.043500e+01 -6.435000e+00 1.443500e+01 1.565000e+00 3.043500e+01 3.356500e+01 9.443500e+01 1.615650e+02 3.504350e+02 6.235440e+02 1.474477e+03 2.325409e+03 6.162705e+03 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04] q1 = -38967.372 W1 = W0 Branch 2 m6 = -0.00905004557884357 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 19 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 2, 4, 5, 7, 8, 9, 11, 3, 13] null space used norm(d0) = 13.73 alpha tilde = 1.06591256 x1 = [7.080000e-01 1.583000e+00 2.417000e+00 5.583000e+00 1.041700e+01 2.158300e+01 3.112600e+01 1.081650e+02 1.478350e+02 3.641650e+02 6.178630e+02 1.472110e+03 2.326356e+03 6.163178e+03 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04] q1 = -38967.423 W1 = W0 Branch 2 m8 = -0.0037368474923263548 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 20 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 2, 4, 5, 7, 9, 11, 3, 13] schur complement used norm(d0) = 15.423 alpha tilde = 1.42652892 x1 = [3.64000e-01 2.27200e+00 1.72800e+00 6.27200e+00 9.72800e+00 2.22720e+01 3.83520e+01 9.30230e+01 1.47694e+02 3.79588e+02 6.11481e+02 1.46945e+03 2.32742e+03 6.16371e+03 1.00000e+04 1.00000e+04 1.00000e+04 1.00000e+04 1.00000e+04 1.00000e+04] q1 = -38967.452 W1 = W0 Branch 2 m4 = -0.0003537463373766392 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 21 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 2, 5, 7, 9, 11, 3, 13] schur complement used norm(d0) = 1.317 alpha tilde = 2.12154936 x1 = [6.640000e-01 1.672000e+00 2.328000e+00 5.672000e+00 9.370000e+00 2.358800e+01 3.780700e+01 9.279800e+01 1.477880e+02 3.796260e+02 6.114650e+02 1.469444e+03 2.327423e+03 6.163711e+03 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04] q1 = -38967.452 W1 = W0 Branch 2 m2 = -3.5290783356346474e-05 < 0. Back to Branch 1. ----------------------------------------------------------- ITERATION 22 ----------------------------------------------------------- Branch 1; k = 0 W0 = [1, 5, 7, 9, 11, 3, 13] schur complement used norm(d0) = 0.145 alpha tilde = 4.73835548 x1 = [7.350000e-01 1.530000e+00 2.326000e+00 5.818000e+00 9.309000e+00 2.356300e+01 3.781700e+01 9.280200e+01 1.477860e+02 3.796260e+02 6.114650e+02 1.469444e+03 2.327422e+03 6.163711e+03 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04] q1 = -38967.452 W1 = W0 Branch 2 mj >= 0 forall j. Solution found. ----------------------------------------------------------- ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND. ----------------------------------------------------------- x* = [7.350000e-01 1.530000e+00 2.326000e+00 5.818000e+00 9.309000e+00 2.356300e+01 3.781700e+01 9.280200e+01 1.477860e+02 3.796260e+02 6.114650e+02 1.469444e+03 2.327422e+03 6.163711e+03 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04] W* = [1, 5, 7, 9, 11, 3, 13]
Resultados¶
print("Optimum reached at x* =",np.round(xf,decimals=3))
print("Function at x* is q(x*) =",np.round(problem_km.evaluate(xf),decimals=3))
print("Final active set is W* =",Wf)
Optimum reached at x* = [7.350000e-01 1.530000e+00 2.326000e+00 5.818000e+00 9.309000e+00 2.356300e+01 3.781700e+01 9.280200e+01 1.477860e+02 3.796260e+02 6.114650e+02 1.469444e+03 2.327422e+03 6.163711e+03 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04] Function at x* is q(x*) = -38967.452 Final active set is W* = [1, 5, 7, 9, 11, 3, 13]
n = q_traj.size
x = np.arange(n)
fig, axs = plt.subplots(1, 1, figsize=(12, 5))
fig.suptitle('q optimization per iteration', fontsize=14)
axs.plot(x, q_traj, marker='o', color='black')
axs.set_xlabel('Step')
axs.set_xticks(x)
axs.set_ylabel('q')
axs.grid(True)
n = int(problem_km.G.shape[1])
x = cp.Variable(n)
prob = cp.Problem(
cp.Minimize((1 / 2) * cp.quad_form(x, problem_km.G) +
problem_km.c.T @ x), [problem_km.A[problem_km.restrictions_l,:] @ x <= problem_km.b[problem_km.restrictions_l]]
)
prob.solve()
solution_cvxpy = x.value
print("Method used:",prob.solver_stats.solver_name)
print("Our optimum value:", round(float(problem_km.evaluate(xf)),2))
print("Optimum value obtained with cvxpy:",round(prob.value,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-solution_cvxpy)/(np.linalg.norm(xf)+np.linalg.norm(solution_cvxpy)),2))
Method used: OSQP Our optimum value: -38967.45 Optimum value obtained with cvxpy: -38967.45 Normalized euclidean distance: 0.0
Observaciones¶
- La solución obtenida con nuestro algoritmo es identica a la obtenida con CVXPY (distancia normalizada cerca de 0).
- Utilizamos x0 inicial factible e imponiendo a W0 como restricciones activas.
print("W0:",W0)
print("iteraciones:",x_traj.size)
print("Valor óptimo obtenido:",round(float(problem_km.evaluate(xf)),2))
W0: [18, 12, 17, 14, 13, 16, 19, 15] iteraciones: 1340 Valor óptimo obtenido: -38967.45
V. Problema 3 - lp_blend¶
A partir de los datos del archivo lp_blend.mat, se plantea el siguiente problema de optimización cuadrática. Dado $G \in \mathbb{R}^{n \times n}$, se desea resolver:
Este planteamiento puede implementarse dentro de una función separada. Una vez definido el problema, se debe aplicar el algoritmo desarrollado para resolverlo, siguiendo los pasos:
- Encontrar un punto inicial $\vec{x_0}$ mediante programación lineal (utilizando
from scipy.optimize import linprog), considerando solo $A$ y $\vec{b}$. - Definir el conjunto activo inicial $W_0$ y un vector $z \in E$.
Se deben documentar los siguientes elementos:
- El conjunto activo inicial $W_0$,
- El número de iteraciones realizadas,
- El valor óptimo alcanzado.
problem_load = loadProblem("lp_blend.mat")
Ae = problem_load["AE"]
n = Ae.shape[1]
Al = -np.identity(n)
A = np.vstack((Ae, Al))
be = problem_load["bE"]
bl = np.zeros(n)
b = np.concatenate((be, bl))
c = problem_load["c"]
G = np.identity(n)
restrictions_eps = np.arange(0, Ae.shape[0]).tolist()
restrictions_l = np.arange(Ae.shape[0], Ae.shape[0] + n).tolist()
problem_blend = QuadraticProgramming(G, c, A, b, restrictions_eps, restrictions_l)
x0 = problem_blend.generate_feasible()
W0 = restrictions_eps.copy()
res = active_set(problem_blend, x0, W0, trajectory=True, tol=1e-9)
xf, wf, x_traj, q_traj = res
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf] ----------------------------------------------------------- ITERATION 0 ----------------------------------------------------------- Branch 1; k = 0 W0 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73] null space used norm(d0) = 14.405 alpha tilde = 1.47629031 x1 = [-3.6200e-01 -1.3760e+00 1.6860e+00 2.3000e-02 1.5700e-01 -4.8600e-01 5.4000e-02 -3.5400e-01 2.0320e+00 5.1000e-02 1.2100e-01 6.8000e-02 2.6000e-01 7.3000e-02 1.9570e+00 6.0000e-02 -3.1000e-02 1.8200e-01 7.6500e-01 1.7830e+00 1.6510e+00 8.8550e+00 1.6940e+00 1.5243e+01 1.2154e+01 1.0247e+01 3.5510e+00 4.6240e+00 1.1863e+01 1.4000e-01 -9.8000e-02 1.1077e+01 8.8960e+00 5.9480e+00 2.6070e+00 -3.2500e-01 2.5000e-02 1.1700e-01 3.1750e+00 5.9000e-02 -6.2800e-01 1.9310e+00 1.2440e+00 1.9900e+00 6.1600e-01 -6.9500e-01 -3.0300e-01 -1.4420e+00 1.4690e+00 1.9180e+00 2.3090e+00 7.7600e-01 3.7300e-01 -1.2360e+00 -6.2700e-01 5.9700e-01 6.7300e-01 4.1400e-01 -1.2000e-01 3.6500e-01 1.1840e+00 4.6300e-01 1.1810e+00 -1.2360e+00 -6.2700e-01 1.9000e-02 1.9990e+00 1.2280e+00 4.3700e-01 -2.7100e-01 3.2500e-01 4.1000e-02 7.4700e-01 6.1000e-02 1.9890e+00 1.3400e+00 8.0500e-01 -4.2500e-01 1.1200e-01 2.4000e-02 7.0400e-01 1.2000e-02 1.7380e+00 1.2320e+00 1.7300e-01 1.0450e+00 6.4100e-01 1.6850e+00 4.4700e-01 2.9870e+00 3.4350e+00 2.5230e+00 2.1190e+00 4.6420e+00 4.1000e-01 1.6010e+00 1.5500e+00 2.7800e-01 4.1000e-02 1.1700e-01 3.9960e+00 4.2050e+00 -1.3400e-01 -5.5300e-01 -1.3400e-01 4.6300e-01 1.1600e-01 -2.5000e-01 2.5200e-01 -2.7830e+00 -4.3000e-01 3.0000e-02 4.6200e+00 1.0070e+00] q1 = 587.078 W1 = W0 Branch 2 mj >= 0 forall j. Solution found. ----------------------------------------------------------- ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND. ----------------------------------------------------------- x* = [-3.6200e-01 -1.3760e+00 1.6860e+00 2.3000e-02 1.5700e-01 -4.8600e-01 5.4000e-02 -3.5400e-01 2.0320e+00 5.1000e-02 1.2100e-01 6.8000e-02 2.6000e-01 7.3000e-02 1.9570e+00 6.0000e-02 -3.1000e-02 1.8200e-01 7.6500e-01 1.7830e+00 1.6510e+00 8.8550e+00 1.6940e+00 1.5243e+01 1.2154e+01 1.0247e+01 3.5510e+00 4.6240e+00 1.1863e+01 1.4000e-01 -9.8000e-02 1.1077e+01 8.8960e+00 5.9480e+00 2.6070e+00 -3.2500e-01 2.5000e-02 1.1700e-01 3.1750e+00 5.9000e-02 -6.2800e-01 1.9310e+00 1.2440e+00 1.9900e+00 6.1600e-01 -6.9500e-01 -3.0300e-01 -1.4420e+00 1.4690e+00 1.9180e+00 2.3090e+00 7.7600e-01 3.7300e-01 -1.2360e+00 -6.2700e-01 5.9700e-01 6.7300e-01 4.1400e-01 -1.2000e-01 3.6500e-01 1.1840e+00 4.6300e-01 1.1810e+00 -1.2360e+00 -6.2700e-01 1.9000e-02 1.9990e+00 1.2280e+00 4.3700e-01 -2.7100e-01 3.2500e-01 4.1000e-02 7.4700e-01 6.1000e-02 1.9890e+00 1.3400e+00 8.0500e-01 -4.2500e-01 1.1200e-01 2.4000e-02 7.0400e-01 1.2000e-02 1.7380e+00 1.2320e+00 1.7300e-01 1.0450e+00 6.4100e-01 1.6850e+00 4.4700e-01 2.9870e+00 3.4350e+00 2.5230e+00 2.1190e+00 4.6420e+00 4.1000e-01 1.6010e+00 1.5500e+00 2.7800e-01 4.1000e-02 1.1700e-01 3.9960e+00 4.2050e+00 -1.3400e-01 -5.5300e-01 -1.3400e-01 4.6300e-01 1.1600e-01 -2.5000e-01 2.5200e-01 -2.7830e+00 -4.3000e-01 3.0000e-02 4.6200e+00 1.0070e+00] W* = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]
Resultados¶
print("Optimum reached at x* =",np.round(xf,decimals=3))
print("Function at x* is q(x*) =",np.round(problem_blend.evaluate(xf),decimals=3))
print("Final active set is W* =",Wf)
Optimum reached at x* = [-3.6200e-01 -1.3760e+00 1.6860e+00 2.3000e-02 1.5700e-01 -4.8600e-01 5.4000e-02 -3.5400e-01 2.0320e+00 5.1000e-02 1.2100e-01 6.8000e-02 2.6000e-01 7.3000e-02 1.9570e+00 6.0000e-02 -3.1000e-02 1.8200e-01 7.6500e-01 1.7830e+00 1.6510e+00 8.8550e+00 1.6940e+00 1.5243e+01 1.2154e+01 1.0247e+01 3.5510e+00 4.6240e+00 1.1863e+01 1.4000e-01 -9.8000e-02 1.1077e+01 8.8960e+00 5.9480e+00 2.6070e+00 -3.2500e-01 2.5000e-02 1.1700e-01 3.1750e+00 5.9000e-02 -6.2800e-01 1.9310e+00 1.2440e+00 1.9900e+00 6.1600e-01 -6.9500e-01 -3.0300e-01 -1.4420e+00 1.4690e+00 1.9180e+00 2.3090e+00 7.7600e-01 3.7300e-01 -1.2360e+00 -6.2700e-01 5.9700e-01 6.7300e-01 4.1400e-01 -1.2000e-01 3.6500e-01 1.1840e+00 4.6300e-01 1.1810e+00 -1.2360e+00 -6.2700e-01 1.9000e-02 1.9990e+00 1.2280e+00 4.3700e-01 -2.7100e-01 3.2500e-01 4.1000e-02 7.4700e-01 6.1000e-02 1.9890e+00 1.3400e+00 8.0500e-01 -4.2500e-01 1.1200e-01 2.4000e-02 7.0400e-01 1.2000e-02 1.7380e+00 1.2320e+00 1.7300e-01 1.0450e+00 6.4100e-01 1.6850e+00 4.4700e-01 2.9870e+00 3.4350e+00 2.5230e+00 2.1190e+00 4.6420e+00 4.1000e-01 1.6010e+00 1.5500e+00 2.7800e-01 4.1000e-02 1.1700e-01 3.9960e+00 4.2050e+00 -1.3400e-01 -5.5300e-01 -1.3400e-01 4.6300e-01 1.1600e-01 -2.5000e-01 2.5200e-01 -2.7830e+00 -4.3000e-01 3.0000e-02 4.6200e+00 1.0070e+00] Function at x* is q(x*) = 587.078 Final active set is W* = [1, 5, 7, 9, 11, 3, 13]
n = q_traj.size
x = np.arange(n)
plt.plot(x, q_traj, marker='o',color='black')
plt.xlabel('Step')
plt.xticks(x)
plt.ylabel('q')
plt.title('q optimization per iteration')
plt.grid(True)
plt.show()
Comparamos solución con CVXPY.
# Define and solve the CVXPY problem.
problem_load = loadProblem("lp_blend.mat")
Ae = problem_load["AE"]
n = Ae.shape[1]
Al = -np.identity(n)
A = np.vstack((Ae, Al))
be = problem_load["bE"]
bl = np.zeros(n)
b = np.concatenate((be, bl))
c = problem_load["c"]
G = np.identity(n)
x = cp.Variable(n)
prob = cp.Problem(
cp.Minimize((1 / 2) * cp.quad_form(x, G) +
c.T @ x), [Al @ x <= bl, Ae @ x == be]
)
prob.solve()
solution = x.value
print("Method used:",prob.solver_stats.solver_name)
print("Active set optimum value:", round(float(problem_blend.evaluate(xf)),2))
print("CVXPY optimum value:",round(prob.value,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-x.value)/(np.linalg.norm(xf)+np.linalg.norm(x.value)),2))
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf] Method used: OSQP Active set optimum value: 587.08 CVXPY optimum value: 712.13 Normalized euclidean distance: 0.44
Nótese que CVXPY utiliza el método OSQP. Los valores resultados obtenidos con ambos métodos son similares (distancia normalizada < 0.5). El método del conjunto activo nos da un mejor minimizador.
Observaciones adicionales¶
Al trabajar con matrices dispersas en optimización con igualdades, pueden surgir problemas de precisión numérica, afectando la estabilidad de la solución, especialmente si las matrices son mal condicionadas o las tolerancias no son adecuadas, lo que puede generar resultados ligeramente diferentes dependiendo del solver y la configuración numérica.
print("x0:",x0)
print("W0:",W0)
print("iteraciones:",x_traj.size)
print("Valor óptimo obtenido:",round(float(problem_blend.evaluate(xf)),2))
x0: [-0. -0. 0. 0. -0. -0. -0. -0. 0. 0. -0. 0. -0. -0. 0. 0. -0. -0. 0. -0. 0. 23.26 5.25 26.32 21.05 13.45 2.58 10. 10. 0. 0. 0. -0. -0. -0. -0. -0. -0. -0. -0. 0. -0. -0. -0. -0. 0. 0. 0. -0. -0. -0. 0. 0. -0. -0. -0. -0. 0. 0. -0. 0. -0. 0. 0. -0. -0. -0. -0. 0. 0. -0. 0. -0. -0. -0. -0. 0. 0. 0. 0. -0. -0. -0. -0. -0. -0. 0. -0. -0. -0. -0. -0. -0. -0. -0. -0. 0. -0. -0. -0. -0. 0. 0. 0. 0. 0. -0. 0. -0. 0. 0. -0. -0. -0. ] W0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73] iteraciones: 228 Valor óptimo obtenido: 587.08
VI. Referencias¶
[1] Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. Springer.
[2] Andreas Wachtel, Curso de Optimización Numérica, Instituto Tecnológico Autónomo de México (ITAM), primavera 2025.
[3] Harris, C. R., Millman, K. J., van der Walt, S. J., et al. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. https://doi.org/10.1038/s41586-020-2649-2
[4] Virtanen, P., Gommers, R., Oliphant, T. E., et al. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2
[5] Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9(3), 90–95. https://doi.org/10.1109/MCSE.2007.55
[6] Caswell, T. A., Droettboom, M., Lee, A., et al. (2016). Matplotlib/matplotlib: v1.5.1. Zenodo. https://doi.org/10.5281/zenodo.44579
[7] Diamond, S., & Boyd, S. (2016). CVXPY: A Python-Embedded Modeling Language for Convex Optimization. Journal of Machine Learning Research, 17(83), 1–5. http://jmlr.org/papers/v17/14-541.html
[8] Plotly Technologies Inc. (2015). Collaborative data science. Montreal, QC. https://plotly.com/python/
[9] Python Software Foundation. random — Generate pseudo-random numbers. Documentación oficial de Python 3. https://docs.python.org/3/library/random.html
[10] Andreas Wachtel, Instituto Tecnológico Autónomo de México (ITAM), primavera 2025. Archivo AW_loadProblem.py.